/*
 * =====================================================================================
 *
 *       Filename:  skeletonization.cpp
 *
 *    Description:  Displays skeletonization of character
 *
 *        Version:  1.0
 *        Created:  05/07/12 15:33:25
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Moos Hueting (mh), mooshueting@gmail.com
 *   Organization:  
 *
 * =====================================================================================
 */
#include <iostream>
#include <cv.h>
#include <highgui.h>

using namespace std;
using namespace cv;

Mat Skeletonization(Mat *src) {
    Mat skel(src->size(), CV_8UC1, Scalar(0));
    Mat temp(src->size(), CV_8UC1);

    Mat element = getStructuringElement(MORPH_CROSS, Size(3, 3));

    bool done;
    do {
        morphologyEx(*src, temp, cv::MORPH_OPEN, element);
        bitwise_not(temp, temp);
        bitwise_and(*src, temp, temp);
        bitwise_or(skel, temp, skel);
        erode(*src, *src, element);

        done = (cv::norm(*src) == 0);
    } while (!(done));

    return skel;

}

void ThinSubiteration1(CvMat *pSrc, CvMat *pDst) {
	int rows = pSrc->rows;
	int cols = pSrc->cols;
	cvCopy(pSrc, pDst);
	for(int i = 0; i < rows; i++) {
		for(int j = 0; j < cols; j++) {
			if(CV_MAT_ELEM(*pSrc, float, i, j) == 1.0f) {
				/// get 8 neighbors
				/// calculate C(p)
				int neighbor0 = (int) CV_MAT_ELEM(*pSrc, float, i-1, j-1);
				int neighbor1 = (int) CV_MAT_ELEM(*pSrc, float, i-1, j);
				int neighbor2 = (int) CV_MAT_ELEM(*pSrc, float, i-1, j+1);
				int neighbor3 = (int) CV_MAT_ELEM(*pSrc, float, i, j+1);
				int neighbor4 = (int) CV_MAT_ELEM(*pSrc, float, i+1, j+1);
				int neighbor5 = (int) CV_MAT_ELEM(*pSrc, float, i+1, j);
				int neighbor6 = (int) CV_MAT_ELEM(*pSrc, float, i+1, j-1);
				int neighbor7 = (int) CV_MAT_ELEM(*pSrc, float, i, j-1);
				int C = int(~neighbor1 & ( neighbor2 | neighbor3)) +
						 int(~neighbor3 & ( neighbor4 | neighbor5)) +
						 int(~neighbor5 & ( neighbor6 | neighbor7)) +
						 int(~neighbor7 & ( neighbor0 | neighbor1));
				if(C == 1) {
					/// calculate N
					int N1 = int(neighbor0 | neighbor1) +
							 int(neighbor2 | neighbor3) +
							 int(neighbor4 | neighbor5) +
							 int(neighbor6 | neighbor7);
					int N2 = int(neighbor1 | neighbor2) +
							 int(neighbor3 | neighbor4) +
							 int(neighbor5 | neighbor6) +
							 int(neighbor7 | neighbor0);
					int N = min(N1,N2);
					if ((N == 2) || (N == 3)) {
						/// calculate criteria 3
						int c3 = ( neighbor1 | neighbor2 | ~neighbor4) & neighbor3;	
						if(c3 == 0) {
							CV_MAT_ELEM(*pDst, float, i, j) = 0.0f;
						}
					}
				}
			}
		}
	}
}

void ThinSubiteration2(CvMat *pSrc, CvMat *pDst) {
	int rows = pSrc->rows;
	int cols = pSrc->cols;
	cvCopy(pSrc, pDst);
	for(int i = 0; i < rows; i++) {
		for(int j = 0; j < cols; j++) {
			if ( CV_MAT_ELEM(*pSrc, float, i, j) == 1.0f) {
				/// get 8 neighbors
				/// calculate C(p)
				int neighbor0 = (int) CV_MAT_ELEM(*pSrc, float, i-1, j-1);
				int neighbor1 = (int) CV_MAT_ELEM(*pSrc, float, i-1, j);
				int neighbor2 = (int) CV_MAT_ELEM(*pSrc, float, i-1, j+1);
				int neighbor3 = (int) CV_MAT_ELEM(*pSrc, float, i, j+1);
				int neighbor4 = (int) CV_MAT_ELEM(*pSrc, float, i+1, j+1);
				int neighbor5 = (int) CV_MAT_ELEM(*pSrc, float, i+1, j);
				int neighbor6 = (int) CV_MAT_ELEM(*pSrc, float, i+1, j-1);
				int neighbor7 = (int) CV_MAT_ELEM(*pSrc, float, i, j-1);
				int C = int(~neighbor1 & ( neighbor2 | neighbor3)) +
					int(~neighbor3 & ( neighbor4 | neighbor5)) +
					int(~neighbor5 & ( neighbor6 | neighbor7)) +
					int(~neighbor7 & ( neighbor0 | neighbor1));
				if(C == 1) {
					/// calculate N
					int N1 = int(neighbor0 | neighbor1) +
						int(neighbor2 | neighbor3) +
						int(neighbor4 | neighbor5) +
						int(neighbor6 | neighbor7);
					int N2 = int(neighbor1 | neighbor2) +
						int(neighbor3 | neighbor4) +
						int(neighbor5 | neighbor6) +
						int(neighbor7 | neighbor0);
					int N = min(N1,N2);
					if((N == 2) || (N == 3)) {
						int E = (neighbor5 | neighbor6 | ~neighbor0) & neighbor7;
						if(E == 0) {
							CV_MAT_ELEM(*pDst, float, i, j) = 0.0f;
						}
					}
				}
			}
		}
	}
}

void MorphologicalThinning(CvMat *pSrc, CvMat *pDst) {
	bool bDone = false;
	int rows = pSrc->rows;
	int cols = pSrc->cols;
	/// pad source
	CvMat *p_enlarged_src = cvCreateMat(rows + 2, cols + 2, CV_32FC1);
	for(int i = 0; i < (rows+2); i++) {
		CV_MAT_ELEM(*p_enlarged_src, float, i, 0)		= 0.0f;
		CV_MAT_ELEM(*p_enlarged_src, float, i, cols+1)	= 0.0f;
	}
	for(int j = 0; j < (cols+2); j++) {
		CV_MAT_ELEM(*p_enlarged_src, float, 0, j)		= 0.0f;
		CV_MAT_ELEM(*p_enlarged_src, float, rows+1, j)	= 0.0f;
	}
	for(int i = 0; i < rows; i++) {
		for(int j = 0; j < cols; j++) {
			if (CV_MAT_ELEM(*pSrc, float, i, j) >= 0.5f) {
				CV_MAT_ELEM(*p_enlarged_src, float, i+1, j+1) = 1.0f;
			}
			else
				CV_MAT_ELEM(*p_enlarged_src, float, i+1, j+1) = 0.0f;
		}
	}
	/// start to thin
	CvMat *p_thinMat1	= cvCreateMat(rows + 2, cols + 2, CV_32FC1);
	CvMat *p_thinMat2	= cvCreateMat(rows + 2, cols + 2, CV_32FC1);
	CvMat *p_cmp		= cvCreateMat(rows + 2, cols + 2, CV_8UC1);
	while (bDone != true) {
		/// sub-iteration 1
		ThinSubiteration1(p_enlarged_src, p_thinMat1);
		/// sub-iteration 2
		ThinSubiteration2(p_thinMat1, p_thinMat2);
		/// compare
		cvCmp(p_enlarged_src, p_thinMat2, p_cmp, CV_CMP_EQ);
		/// check
		int num_non_zero = cvCountNonZero(p_cmp);
		if(num_non_zero == (rows + 2) * (cols + 2)) {
			bDone = true;
		}
		/// copy
		cvCopy(p_thinMat2, p_enlarged_src);
	}
	/// copy result
	for(int i = 0; i < rows; i++) {
		for(int j = 0; j < cols; j++) {
			CV_MAT_ELEM(*pDst, float, i, j) = CV_MAT_ELEM(*p_enlarged_src, float, i+1, j+1);
		}
	}
	/// clean memory
	cvReleaseMat(&p_enlarged_src);
	cvReleaseMat(&p_thinMat1);
	cvReleaseMat(&p_thinMat2);
	cvReleaseMat(&p_cmp);
}

int main(int argc, char* argv[]) {

    Mat image, imageResized;

    image = imread(argv[1], CV_LOAD_IMAGE_GRAYSCALE );
    threshold(image, image, 100, 255, THRESH_BINARY );

    resize(image, imageResized, Size(0,0), 10, 10, INTER_NEAREST);

    namedWindow("Display", CV_WINDOW_AUTOSIZE);
    imshow("Display", imageResized);
    waitKey(0);

    image = Skeletonization(&image);

    resize(image, imageResized, Size(0,0), 10, 10, INTER_NEAREST);

    imshow("Display", imageResized);
    waitKey(0);

}
